{
"cells": [
{
"cell_type": "markdown",
"id": "cell-0",
"metadata": {},
"source": [
"# SGP4 vs Numerical Propagation\n",
"\n",
"SGP4 and numerical propagation are two fundamentally different approaches to predicting satellite orbits:\n",
"\n",
"- **SGP4** is an analytical model that uses simplified perturbation theory (J2--J4 zonal harmonics + atmospheric drag via B*). It is fast and standardized, but limited in accuracy — TLEs are designed to be used *only* with SGP4, and their accuracy degrades over days.\n",
"- **Numerical propagation** integrates the full equations of motion with configurable force models (high-order gravity, lunisolar perturbations, drag with NRLMSISE-00, solar radiation pressure). It is slower but can achieve meter-level accuracy.\n",
"\n",
"This tutorial propagates the same orbit both ways and shows how and why they diverge."
]
},
{
"cell_type": "markdown",
"id": "cell-1",
"metadata": {},
"source": [
"## Setup: Common Initial Conditions\n",
"\n",
"We start from an ISS TLE. SGP4 propagates the TLE directly, outputting position in the TEME frame. For the numerical propagator, we convert the SGP4 state at TLE epoch from TEME to GCRF."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-2",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"\n",
"# ISS TLE\n",
"tle = sk.TLE.from_lines([\n",
" \"ISS (ZARYA)\",\n",
" \"1 25544U 98067A 24100.50000000 .00016717 00000-0 30233-3 0 9993\",\n",
" \"2 25544 51.6416 247.4627 0006703 130.5360 229.6208 15.49438695449041\"\n",
"])\n",
"\n",
"epoch = tle.epoch\n",
"print(f\"TLE epoch: {epoch}\")\n",
"print(f\"Mean motion: {tle.mean_motion:.8f} rev/day\")\n",
"print(f\"Inclination: {tle.inclination:.4f} deg\")\n",
"print(f\"B*: {tle.bstar:.6e}\")\n",
"\n",
"# SGP4 state at epoch (TEME frame)\n",
"pos_teme, vel_teme = sk.sgp4(tle, epoch)\n",
"print(\"\\nSGP4 at epoch (TEME):\")\n",
"print(f\" pos = [{pos_teme[0]/1e3:.3f}, {pos_teme[1]/1e3:.3f}, {pos_teme[2]/1e3:.3f}] km\")\n",
"print(f\" vel = [{vel_teme[0]/1e3:.6f}, {vel_teme[1]/1e3:.6f}, {vel_teme[2]/1e3:.6f}] km/s\")\n",
"print(f\" altitude ≈ {(np.linalg.norm(pos_teme) - sk.consts.earth_radius) / 1e3:.1f} km\")\n",
"\n",
"# Convert to GCRF for numerical propagator\n",
"q_teme2gcrf = sk.frametransform.qteme2gcrf(epoch)\n",
"pos_gcrf = q_teme2gcrf * pos_teme\n",
"vel_gcrf = q_teme2gcrf * vel_teme\n",
"state0 = np.concatenate([pos_gcrf, vel_gcrf])\n",
"print(\"\\nInitial state (GCRF):\")\n",
"print(f\" pos = [{pos_gcrf[0]/1e3:.3f}, {pos_gcrf[1]/1e3:.3f}, {pos_gcrf[2]/1e3:.3f}] km\")"
]
},
{
"cell_type": "markdown",
"id": "cell-3",
"metadata": {},
"source": [
"## Propagate Both for 3 Days\n",
"\n",
"We propagate both models forward and compare positions in the GCRF frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-4",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"\n",
"# Time array: 3 days, sampled every 2 minutes\n",
"n_days = 3\n",
"dt_minutes = 2\n",
"n_points = int(n_days * 24 * 60 / dt_minutes) + 1\n",
"times = [epoch + sk.duration(minutes=i * dt_minutes) for i in range(n_points)]\n",
"\n",
"# --- SGP4 propagation ---\n",
"sgp4_teme_pos, sgp4_teme_vel = sk.sgp4(tle, times)\n",
"\n",
"# Convert all SGP4 results from TEME to GCRF\n",
"q_arr = [sk.frametransform.qteme2gcrf(t) for t in times]\n",
"sgp4_gcrf = np.array([q * p for q, p in zip(q_arr, sgp4_teme_pos)])\n",
"\n",
"# --- Numerical propagation ---\n",
"# High-fidelity force model\n",
"settings = sk.propsettings(\n",
" gravity_model=sk.gravmodel.jgm3,\n",
" gravity_degree=20,\n",
" abs_error=1e-10,\n",
" rel_error=1e-10,\n",
")\n",
"\n",
"result = sk.propagate(\n",
" state0,\n",
" epoch,\n",
" end=epoch + sk.duration(days=n_days),\n",
" propsettings=settings,\n",
")\n",
"\n",
"# Interpolate numerical solution at the same times\n",
"num_gcrf = np.array([result.interp(t)[:3] for t in times])\n",
"\n",
"# --- Position difference ---\n",
"pos_diff = np.linalg.norm(sgp4_gcrf - num_gcrf, axis=1)\n",
"hours = np.array([(t - epoch).hours for t in times])\n",
"\n",
"print(\"Position difference:\")\n",
"print(f\" After 1 hour: {pos_diff[30]:.1f} m\")\n",
"print(f\" After 6 hours: {pos_diff[180]:.1f} m\")\n",
"print(f\" After 24 hours: {pos_diff[720]:.0f} m\")\n",
"print(f\" After 72 hours: {pos_diff[-1]:.0f} m\")"
]
},
{
"cell_type": "markdown",
"id": "cell-5",
"metadata": {},
"source": [
"## Visualize the Divergence"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-6",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']\n",
"\n",
"fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)\n",
"\n",
"ax = axes[0]\n",
"ax.plot(hours, pos_diff / 1e3, linewidth=1.5)\n",
"ax.set_ylabel(\"Position Difference [km]\")\n",
"ax.set_title(\"SGP4 vs Numerical Propagation (ISS)\")\n",
"\n",
"# Component-wise differences\n",
"ax = axes[1]\n",
"diff = sgp4_gcrf - num_gcrf\n",
"ax.plot(hours, diff[:, 0] / 1e3, label=\"X\", linewidth=0.8)\n",
"ax.plot(hours, diff[:, 1] / 1e3, label=\"Y\", linewidth=0.8)\n",
"ax.plot(hours, diff[:, 2] / 1e3, label=\"Z\", linewidth=0.8)\n",
"ax.set_xlabel(\"Time Since Epoch [hours]\")\n",
"ax.set_ylabel(\"Component Difference [km]\")\n",
"ax.set_title(\"GCRF Component Differences\")\n",
"ax.legend()\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "cell-7",
"metadata": {},
"source": [
"## Effect of Gravity Model Fidelity\n",
"\n",
"SGP4 only models J2--J4 (zonal harmonics up to degree 4). The numerical propagator can use much higher degree/order. Let's compare the divergence with different gravity truncation levels to see how much of the error comes from higher-order gravity terms."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-8",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"\n",
"gravity_degrees = [4, 10, 20, 40]\n",
"results_by_degree = {}\n",
"\n",
"for deg in gravity_degrees:\n",
" settings = sk.propsettings(\n",
" gravity_model=sk.gravmodel.jgm3,\n",
" gravity_degree=deg,\n",
" abs_error=1e-10,\n",
" rel_error=1e-10,\n",
" )\n",
" res = sk.propagate(state0, epoch,\n",
" end=epoch + sk.duration(days=n_days),\n",
" propsettings=settings)\n",
" num_pos = np.array([res.interp(t)[:3] for t in times])\n",
" results_by_degree[deg] = np.linalg.norm(sgp4_gcrf - num_pos, axis=1)\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 5))\n",
"for deg in gravity_degrees:\n",
" ax.plot(hours, results_by_degree[deg] / 1e3, linewidth=1.2, label=f\"Degree {deg}\")\n",
"ax.set_xlabel(\"Time Since Epoch [hours]\")\n",
"ax.set_ylabel(\"Position Difference [km]\")\n",
"ax.set_title(\"SGP4 vs Numerical: Effect of Gravity Model Order\")\n",
"ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "cell-9",
"metadata": {},
"source": [
"## Altitude Comparison\n",
"\n",
"Another way to see the divergence: compare the altitude profiles. SGP4's simplified drag model (B*) produces slightly different orbital decay than the NRLMSISE-00 density model used by the numerical propagator."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-10",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"\n",
"# SGP4 altitudes (via ITRF -> geodetic)\n",
"sgp4_alt = []\n",
"for t, p_teme in zip(times, sgp4_teme_pos):\n",
" p_itrf = sk.frametransform.qteme2itrf(t) * p_teme\n",
" sgp4_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3)\n",
"\n",
"# Numerical altitudes (GCRF -> ITRF -> geodetic)\n",
"num_alt = []\n",
"for t in times:\n",
" state = result.interp(t)\n",
" p_itrf = sk.frametransform.qgcrf2itrf(t) * state[:3]\n",
" num_alt.append(sk.itrfcoord(p_itrf).altitude / 1e3)\n",
"\n",
"fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)\n",
"\n",
"ax = axes[0]\n",
"ax.plot(hours, sgp4_alt, linewidth=0.5, alpha=0.7, label=\"SGP4\")\n",
"ax.plot(hours, num_alt, linewidth=0.5, alpha=0.7, label=\"Numerical\")\n",
"ax.set_ylabel(\"Altitude [km]\")\n",
"ax.set_title(\"Altitude Comparison\")\n",
"ax.legend()\n",
"\n",
"ax = axes[1]\n",
"alt_diff = np.array(sgp4_alt) - np.array(num_alt)\n",
"ax.plot(hours, alt_diff * 1e3, linewidth=0.8)\n",
"ax.set_xlabel(\"Time Since Epoch [hours]\")\n",
"ax.set_ylabel(\"Altitude Difference [m]\")\n",
"ax.set_title(\"SGP4 $-$ Numerical Altitude Difference\")\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}